!
! Copyright (C) 2002-2008 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
MODULE upf_nml

  USE kinds, ONLY: DP
  use radial_grids, ONLY: radial_grid_type

  IMPLICIT NONE
  SAVE

  integer, parameter :: nmesh_max = 2000
  integer, parameter :: nbeta_max = 5
  integer, parameter :: nwfc_max  = 5

     CHARACTER(LEN=80):: generated=' '! generator software
     CHARACTER(LEN=80):: author=' '   ! pseudopotential's author
     CHARACTER(LEN=80):: date=' '     ! generation date
     CHARACTER(LEN=80):: comment=' '  ! author's comment
     CHARACTER(LEN=2) :: psd=' '      ! Element label
     CHARACTER(LEN=20) :: typ=' '     ! Pseudo type ( NC or US or PAW)
     CHARACTER(len=6) :: rel=' '      ! relativistic: {no|scalar|full}
     LOGICAL :: tvanp              ! .true. if Ultrasoft
     LOGICAL :: tcoulombp          ! .true. if Coulomb 1/r potential
     LOGICAL :: nlcc               ! Non linear core corrections
     CHARACTER(LEN=25) :: dft      ! Exch-Corr type
     REAL(DP) :: zp                ! z valence
     REAL(DP) :: etotps            ! total energy
     REAL(DP) :: ecutwfc           ! suggested cut-off for wfc
     REAL(DP) :: ecutrho           ! suggested cut-off for rho
     !
     CHARACTER(len=11) :: nv       ! UPF file three-digit version i.e. 2.0.0
     INTEGER :: lmax               ! maximum l component in beta
     INTEGER :: lmax_rho           ! max l component in charge (should be 2*lmax)

     ! only for single-channel NC PP
     ! Wavefunctions and projectors
     INTEGER :: nwfc               ! number of atomic wavefunctions
     INTEGER :: nbeta              ! number of projectors
     INTEGER :: kbeta( nbeta_max ) ! kbeta(nbeta) see below
     !  kbeta<=mesh is the number of grid points for each beta function
     !              beta(r,nb) = 0 for r > r(kbeta(nb))

     INTEGER :: lll( nbeta_max )     ! lll(nbeta) l of each projector
     REAL(DP) :: beta( nmesh_max, nbeta_max)  ! beta(mesh,nbeta) projectors
     !
     CHARACTER(LEN=2) :: els( nwfc_max )  ! els(nwfc) label of wfc
     CHARACTER(LEN=2) :: els_beta( nbeta_max)  ! els(nbeta) label of beta
     INTEGER  :: nchi( nwfc_max )    ! lchi(nwfc) value of pseudo-n for wavefcts
     INTEGER  :: lchi( nwfc_max )    ! lchi(nwfc) value of l for wavefcts
     REAL(DP) :: oc( nwfc_max )      ! oc(nwfc) occupancies for wavefcts
     REAL(DP) :: epseu( nwfc_max )   ! pseudo one-particle energy (nwfc)
     REAL(DP) :: rcut_chi( nwfc_max )! rcut_chi(nwfc) cutoff inner radius
     REAL(DP) :: rcutus_chi( nwfc_max )! rcutus_chi(nwfc) ultrasoft outerradius
     ! Chi and rho_at are only used for initial density and initial wfcs:
     REAL(DP) :: chi(nmesh_max,nwfc_max)   ! chi(mesh,nwfc) atomic wavefcts
     REAL(DP) :: rho_at(nmesh_max)  ! rho_at(mesh) atomic charge
     ! Minimal radial grid:
     INTEGER :: nmesh               ! number of points in the radial mesh
     REAL(DP) :: xmin              ! the minimum x of the linear mesh
     REAL(DP) :: rmax              ! the maximum radius of the mesh
     REAL(DP) :: zmesh             ! the nuclear charge used for mesh
     REAL(DP) :: dx                ! the deltax of the linear mesh
     REAL(DP) :: r(nmesh_max)     ! r(mesh)  radial grid
     REAL(DP) :: rab(nmesh_max)   ! rab(mesh) dr(x)/dx (x=linear grid)
     ! Pseudized core charge
     REAL(DP) :: rho_atc(nmesh_max) ! rho_atc(mesh) atomic core charge
     ! Local potential
     INTEGER :: lloc                 ! L of channel used to generate local potential
     ! (if < 0 it was generated by smoothing AE potential)
     REAL(DP) :: rcloc               ! vloc = v_ae for r > rcloc
     REAL(DP) :: vloc(nmesh_max)    ! vloc(mesh) local atomic potential
     !
     REAL(DP) :: dion(nbeta_max *nbeta_max) !dion(nbeta,nbeta) atomic D_{mu,nu}

     LOGICAL           :: has_wfc    ! if true, UPF contain AE and PS wfc for each beta

     LOGICAL :: has_so             ! if .true. includes spin-orbit
     REAL(DP) :: rcut(nbeta_max)  ! cut-off radius(nbeta)
     REAL(DP) :: rcutus(nbeta_max)! ultrasoft cut-off radius (nbeta)

     ! PAW:
     LOGICAL  :: tpawp               ! true if atom is PAW, PAW data must be present

     ! GIPAW:
     LOGICAL  :: has_gipaw           ! Whether GIPAW data is included

contains

subroutine upf_nml_read( iunps, upf, grid, ierr )

  !use upf_nml
  use pseudo_types, only: pseudo_upf
  implicit none
  
  integer, intent(in) :: iunps
  type( pseudo_upf ), intent(inout) :: upf
  type( radial_grid_type ), intent(inout), target :: grid
  integer, intent(out) :: ierr

  namelist /header/ &
       generated, author, date, comment, psd, typ, rel, tvanp, tpawp, tcoulombp, &
       has_so, has_wfc, has_gipaw, nlcc, dft, zp, etotps, ecutwfc, ecutrho, nv,  &
       lmax, lmax_rho, lloc, nmesh, nwfc, nbeta

  namelist /mesh/ dx, nmesh, xmin, rmax, zmesh, r, rab

  namelist /nlcc_rho/ rho_atc

  namelist /local/ vloc

  namelist /nonlocal/ els_beta, lll, kbeta, rcut, rcutus, beta, dion

  namelist /pswfc/ els, lchi, oc, nchi, epseu, rcut_chi, rcutus_chi, chi

  namelist /rhoatom/ rho_at

  !open( unit=iunps, file=trim(fname), status='old', delim='APOSTROPHE', iostat=ierr )

  read( iunps, nml=header )

  read( iunps, nml=mesh  )

  if( nlcc )  read( iunps, nml=nlcc_rho )

  read( iunps, nml=local  )

  read( iunps, nml=nonlocal )

  read( iunps, nml=pswfc  )

  read( iunps, nml=rhoatom )

  !close( iunps, status='keep' )

  !copy the upf_nml module to the upf data_structure
  call upf_nml_copy( upf, grid )

  return
end subroutine upf_nml_read

subroutine upf_nml_copy( upf, grid )

  !use upf_nml
  use pseudo_types, only: pseudo_upf
  use radial_grids, only: radial_grid_type, allocate_radial_grid
  implicit none
  
  type( pseudo_upf ), intent(inout) :: upf
  type(radial_grid_type), intent(inout), target :: grid

  integer :: nb, nw

  !nullify elements of the upf that are not used by mini_DFT
  nullify(upf%vnl)
  upf%q_with_l = .false.
  upf%nqf = 0
  upf%nqlc = 0
  upf%qqq_eps = 0.0_dp
  nullify( upf%rinner )
  nullify( upf%qqq )
  nullify( upf%qfunc )
  nullify( upf%qfuncl )
  nullify( upf%qfcoef )
  nullify( upf%aewfc )
  nullify( upf%pswfc )
  nullify( upf%nn )
  nullify( upf%jchi )
  nullify( upf%jjj )
  upf%paw_data_format = 0
  !upf%paw
  upf%paw_as_gipaw = .false.
  upf%gipaw_data_format = 0
  upf%gipaw_ncore_orbitals = 0
  nullify( upf%gipaw_core_orbital_n )
  nullify( upf%gipaw_core_orbital_l )
  nullify( upf%gipaw_core_orbital_el )
  nullify( upf%gipaw_core_orbital )
  nullify( upf%gipaw_vlocal_ae )
  nullify( upf%gipaw_vlocal_ps )
  upf%gipaw_wfs_nchannels = 0
  nullify( upf%gipaw_wfs_el )
  nullify( upf%gipaw_wfs_ll )
  nullify( upf%gipaw_wfs_ae )
  nullify( upf%gipaw_wfs_rcut )
  nullify( upf%gipaw_wfs_rcutus )
  nullify( upf%gipaw_wfs_ps )

  !copy data from the upf_nml module to the upf structure
  upf%nv        = nv
  upf%generated = generated
  upf%author    = author
  upf%date      = date
  upf%comment   = comment
  
  upf%psd = psd
  upf%typ = typ
  upf%rel = rel

  upf%tvanp     = tvanp
  upf%tpawp     = tpawp
  upf%tcoulombp = tcoulombp
  
  upf%has_so       = has_so
  upf%has_wfc      = has_wfc
  upf%has_gipaw    = has_gipaw
  !upf%paw_as_gipaw = paw_as_gipaw

  upf%nlcc     = nlcc
  upf%dft      = dft
  upf%zp       = zp
  upf%etotps   = etotps
  upf%ecutwfc  = ecutwfc
  upf%ecutrho  = ecutrho
  upf%lmax     = lmax
  upf%lmax_rho = lmax_rho
  upf%lloc     = lloc
  upf%mesh     = nmesh
  upf%nwfc     = nwfc
  upf%nbeta    = nbeta

  upf%dx = dx
  upf%mesh = nmesh
  upf%xmin = xmin
  upf%rmax = rmax
  upf%zmesh = zmesh

  CALL allocate_radial_grid(grid, upf%mesh)
  grid%dx    = upf%dx
  grid%mesh  = upf%mesh
  grid%xmin  = upf%xmin
  grid%rmax  = upf%rmax
  grid%zmesh = upf%zmesh
  upf%grid => grid
  upf%r    => upf%grid%r
  upf%rab  => upf%grid%rab

  upf%r(1:upf%mesh) = r(1:upf%mesh)
  upf%rab(1:upf%mesh) = rab(1:upf%mesh)

  ! Reconstruct additional grids
  upf%grid%r2 =  upf%r**2
  upf%grid%sqr = sqrt(upf%r)
  upf%grid%rm1 = upf%r**(-1)
  upf%grid%rm2 = upf%r**(-2)
  upf%grid%rm3 = upf%r**(-3)

  !allocate and initialized core density; initialize to zero if .not.nlcc
  allocate( upf%rho_atc( upf%mesh) )
  upf%rho_atc(1:upf%mesh) = 0._dp
  if( upf%nlcc ) upf%rho_atc(1:upf%mesh) = rho_atc(1:upf%mesh)

  !allocate and initialize local potential
  if( .not. upf%tcoulombp ) then
     allocate( upf%vloc(upf%mesh) )
     upf%vloc(1:upf%mesh) = vloc(1:upf%mesh)
  end if

  !allocate and initialize nonlocal potential
  if( upf%tcoulombp ) upf%nbeta = 0
  IF ( upf%nbeta == 0) then

     !dummy allocation for local-only pp
     upf%nqf = 0
     upf%nqlc= 0
     upf%qqq_eps= -1._dp
     upf%kkbeta = 0  
     ALLOCATE( upf%kbeta(1),         &
               upf%lll(1),           &
               upf%beta(upf%mesh,1), &
               upf%dion(1,1),        &
               upf%rinner(1),        &
               upf%qqq(1,1),         &
               upf%qfunc(upf%mesh,1),&
               upf%qfcoef(1,1,1,1),  &
               upf%rcut(1),          &
               upf%rcutus(1),        &
               upf%els_beta(1) )
     
  ELSE

     !allocate space for non-local part
     ALLOCATE( upf%kbeta(upf%nbeta),          &
               upf%lll(upf%nbeta),            &
               upf%beta(upf%mesh, upf%nbeta), &
               upf%dion(upf%nbeta, upf%nbeta),&
               upf%rcut(upf%nbeta),           &
               upf%rcutus(upf%nbeta),         &
               upf%els_beta(upf%nbeta) )

     !read the projectors
     do nb=1, upf%nbeta

        upf%beta(1:upf%mesh,nb) = beta(1:upf%mesh,nb)
        upf%els_beta(nb) = els_beta(nb)
        upf%lll(nb)      = lll(nb)
        upf%kbeta(nb)    = kbeta(nb)
        upf%rcut(nb)     = rcut(nb)
        upf%rcutus(nb)   = rcutus(nb)

     end do
        
     upf%dion(:,:) = reshape( dion, (/ upf%nbeta, upf%nbeta /) )
     
     upf%kkbeta = maxval( upf%kbeta(1:upf%nbeta) )
     
  END IF

  !allocate and initialize chi functions
  ALLOCATE( upf%chi(upf%mesh,upf%nwfc) )
  ALLOCATE( upf%els(upf%nwfc),        &
            upf%oc(upf%nwfc),         &
            upf%lchi(upf%nwfc),       &
            upf%nchi(upf%nwfc),       &
            upf%rcut_chi(upf%nwfc),   &
            upf%rcutus_chi(upf%nwfc), &
            upf%epseu(upf%nwfc)       )

  DO nw=1, upf%nwfc
     upf%chi(1:upf%mesh, nw) = chi(1:upf%mesh,nw)
     upf%els(nw)   = els(nw)
     upf%lchi(nw)  = lchi(nw)
     upf%oc(nw)    = oc(nw)
     upf%nchi(nw)  = nchi(nw)
     upf%epseu(nw) = epseu(nw)
     upf%rcut_chi(nw)   = rcut_chi(nw)
     upf%rcutus_chi(nw) = rcutus_chi(nw)
  ENDDO

  !allocate and initialize atomic density
  allocate( upf%rho_at( upf%mesh ) )
  upf%rho_at(1:upf%mesh) = rho_at(1:upf%mesh)

end subroutine upf_nml_copy


END MODULE upf_nml

